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We model femiions with an attractive interaction in an optical lattice with a single-band Hubbard 
model away from half-filling with on-site attraction U and nearest neighbor hopping t. Our goal 
is to understand the crossover from BCS superfluidity in the weak attraction limit to the BEG 
of molecules in the strong attraction limit, with particular emphasis on how this crossover in an 
optical lattice dilTers from the much better studied continuum problem. We use a large- A*' theory with 
Sp(2A'^) symmetry to study the fluctuations beyond mean field theory. At T = 0, we calculate across 
the crossover various observables, including chemical potential, gap, ground state energy, speed of 
sound and compressibility. The superfiuid density n^ is found to have non-trivial U/t dependence 
in this lattice system. We show that the transition temperature Tc scales with the energy gap in the 
weak coupling limit but crosses over to a t'^ /U scaling in the BEG limit, where phase fiuctuations 
controlled by Us determine Tc- We also find, quite contrary to our expectations, that in the strong 
coupling limit, the large- Af theory gives qualitatively wrong trends for compressibility. A comparison 
with a simple Hartree shifted BCS theory, which takes into account both pairing and Hartree shifts, 
and correctly recovers the atomic limit and the right qualitative trend for compressibility, reveals 
that the large- A'^ theory on the lattice, although considers a larger number of diagrams, is in fact 
inferior to the simpler Hartree shifted BGS theory. The failure of the large- A^ approach is explained 
by noting (i) the importance of Hartree shift in lattice problems, and (ii) inability of the large- A'^ 
approach to treat particle-particle and particle-hole channels at equal footing at the saddle point 
level. 



PACS numbers: 03.75.Ss, 05.30.Fk 



I. INTRODUCTION 

The problem of the BCS-BEC crossover of strongly 
interacting fermions has been well studied in the contin- 
uum both theoretically [THU] and experimentally [13- 
\TS\ . The system smoothly interpolates between a BCS 
state of loosely bound Cooper pairs to a Bose-Einstein 
condensate of tightly bound diatomic molecules. The 
Leggett-BCS mean field theory gives qualitatively cor- 
rect physics at T = across the crossover and methods 
like functional integral formalism [SHZ], self-consistent 
approximations [S] , large- TV expansions [HI |10! and quan- 
tum Monte Carlo [TTJ [T^ have been used to find quan- 
titative corrections. Experiments have demonstrated the 
condensation of molecules in the BEC limit [T31 E] and 
the superfluidity of the system across resonance has been 
observed [T5] . 

The inclusion of a lattice in the system leads to sev- 
eral qualitative differences with the continuum. One of 
the key features distinguishing a lattice system from the 
continuum is the dependence on interaction strength and 
filling fraction of the superfiuid stiffness of the gas even 
at T = 0. This is in contrast to the continuum case, 
where the T = superfiuid stiffness is fixed by the par- 
ticle mass and density due to Galilean invariance. Con- 
sequently, when phase fiuctuations play a dominant role 
in the loss of phase coherence and the superfiuid stiffness 
sets the scale for transition temperature [TB], the above 
mentioned difference between the lattice and continuum 
becomes explicit. A second difference, which is not en- 
tirely unrelated to the previous point, that arises between 
the continuum and the lattice is regarding the effective 



mass of the bound pairs in the BEC limit. In the con- 
tinuum, the mass of the bound pair in the BEC limit is 
simply twice the mass of the fermions and hence does not 
scale with the coupling strength. In contrast, the effec- 
tive mass of the bosons on the lattice becomes increasing 
larger with the strength of the coupling. This is due to 
the fact that the bosons on the lattice can only move 
around by virtual ionization, and hence the correspond- 
ing hopping matrix element for the bosons, calculated 
within a simple perturbation theory, has an energy de- 
nominator equal to the coupling strength. Consequently, 
the boson mass which is inversely proportional to the 
hopping becomes larger with the strength of coupling. 
Thirdly, on a bipartite lattice there is a particle-hole (p-h) 
transformation that puts additional constraints on ther- 
modynamics. Finally, there is an emergence on a lattice 
of a Charge Density Wave order at half-filling that com- 
petes with the superfiuid (pairing) order. This new order 
arises because at half-filling the lattice Hamiltonian has 
a higher SU(2) symmetry in the spin space that is spon- 
taneously broken. Our primary objective for studying 
attractive fermionic atoms on a lattice is therefore to un- 
derstand how the broken translational invariance affects 
various physical quantities across the crossover. More- 
over, there is a growing interest in performing experi- 
ments with ultra-cold fermionic clouds of both ^°K [T7] 
and ^Li [18] atoms in optical lattices, and although the 
entropy in the current experiments needs to be reduced 
by a factor of 3 or higher in order to reduce the temper- 
ature to below Tc [H], we believe that in future these 
experiments should be able to test the findings of our 
current work. 



The paper is organized as follows: In section II we in- 
troduce the Hamiltonian and discuss the p-h constraints 
on a lattice. In section III we discuss a Hartree shifted 
BCS (HBCS) theory that respects the lattice p-h con- 
straints. In section IV we next outline a diagrammatic 
approach to include the effects of quantum fluctuations 
on top of the HBCS theory and discuss the problem in 
this approach. Next, in section V we develop a large- 
A'' formalism that we use in the rest of the paper. We 
discuss the T = results for different properties of the 
system within the large- A^ approach. We conclude this 
section with a comparison between the HBCS and large- 
N, and show that the former theory gives a more correct 
account of the chemical potential and the compressibility 
at T = 0. In section V we calculate the zero tempera- 
ture superfluid stiffness. In Section VI we outline the 
calculation and results for the critical temperature. We 
conclude in section VIII. 



II. HAMILTONIAN AND PARTICLE-HOLE 
CONSTRAINTS 

In this section we first introduce the Hamiltonian that 
describes two kinds of fermions in a lattice and calcu- 
late the strength of interactions for which a two-particle 
bound state appears (unitarity condition). We next de- 
rive a set of p-h constraints imposed on the thermody- 
namics and develop a simple HBCS theory that respects 
these constraints. 



A. Hamiltonian 

The study of the BCS-BEC crossover in the absence 
of an optical lattice uses the divergence of the scattering 
length near a Feshbach resonance to tune the strength 
of the interactions between the fermions. Although this 
same technique has been applied to fermions in optical 
lattices TS], the Hamiltonian that describes this system 
near resonance is poorly understood. This is due to the 
inherent multi-band nature of the system when the (con- 
tinuum) scattering length between the atoms diverges 
PP] . We do not have a separation of energy scales that 
would allow us to study an effective Hamiltonian in a 
single band. However, as we will show below, the lattice 
strongly modifies the scattering properties of fermions re- 
stricted to the lowest band, to the point that it takes a 
finite amount of on-site interaction to form a (molecular) 
bound state. Thus, a Feshbach resonance is not needed 
to achieve a unitary gas in a lattice. 

The Hamiltonian we will study is the single-band at- 
tractive Hubbard Hamiltonian: 
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. (1) 

Here Cj^- is the fermion annihilation operator at site j, the 
pseudo-spin index a =t, i represent the two-hyperfine 



states, t is the hopping matrix element between adjacent 
sites and the summation indices («, j) represent sums over 
nearest-neighbor sites. The on-site attractive coupling is 
given by — [/ with U > 0, and it is assumed that both the 
hopping t and U are much smaller than the inter-band 
gap. Finally, rii^ = cj^Ci^ is the number operator at site 
i of fermions with spin cr, and /i is the chemical potential. 
For simplicity, we will study homogeneous systems; i.e. 
we neglect the effects of the (typically harmonic) exter- 
nal trapping potential, which can eventually be included 
using a local density approximation. Throughout the 
paper, we have set h — ks — I and we shall use the con- 
vention that all 3-momenta sums are summed over the 
first Brillouin zone and then divided by the total number 
of lattice sites. 

The scattering amplitude between fermions in the lat- 
tice can be obtained by summing up all possible interac- 
tion events of fermions with the dispersion relation ob- 
tained from the kinetic energy in (fTj), et = — 2t[cos(fca;a)-|- 
cos(fcj,a) + cos{kza) — 3] (which we conventionally mea- 
sure from the bottom of the band) , where a is the lattice 
constant. The scattering amplitude can be calculated as 
/ = (TO/47r) r(0, 0) where r(q, w) = U/{1 + t/n(q, cj)) is 
the four-point vertex function for a pair of fermions of 
mass m with center of mass momentum q and H(q, w) 
is the corresponding polarization, which in our case (and 
in the limit T — >► 0) is of the form 
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where the integration is over the Brillouin zone. We can 
now see that the condition for a diverging scattering am- 
plitude (i.e. unitarity) in the lattice is |21) : 
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(3) 



For most experiments, the values of U and t can be more 
or less independently chosen. While U is primarily fixed 
by the magnetic field strength and the latter can be cho- 
sen such that one is always far from a Feshbach reso- 
nance, t can be adjusted by tuning the height of optical 
lattice. Therefore, the single Bloch band picture remains 
valid for the purpose of studying the BCS-BEC crossover 
as depicted by Hamiltonian Q. 



B. Particle-Hole Constraints 

Lattice systems have an additional symmetry stem- 
ming from the possibility of describing the physics in 
terms of either particles or holes; the choice of description 
is usually made in order to simplify the resulting Hamil- 
tonian. In the case of the Hamiltonian ( |l2| ) we can obtain 
an exact relationship between a system with n fermions 
(particles) and one with 2 — n fermions (holes). 

Let us for the moment work in the canonical ensem- 
ble and look for the ground state of the Hamiltonian (fTl) 



with the constraint that the number of particles per site 
is n = n^ + n^. If we now perform the particle-hole trans- 
formation c] ^ = ( — l)*(ii^cr 22 , it can be easily verified 
that the kinetic energy term maintains its form with the 
replacement of the c operators with d operators. On the 
other hand, the on-site interaction term (with the site 
index omitted for clarity) transforms as 
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(4) 
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Given that dld^ + d\d, = 2 — n is fixed in the calcu- 
lation, the terms in the second line of Q are constant 
within the Hilbert space of interest. Thus, the Hamito- 
nian maintains its operational form under the particle- 
hole transformation and the ground state wavefunction 
for a system of n particles is related to the ground state 
wavefunction for a system of 2 — n particles. Their cor- 
responding energies are related as 



£[n) + --=£{2-n) 
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(5) 



Differentiating with respect to n we see that the chemical 
potential, defined as /i(n) = d£{n)/d{n), satisfies 
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(6) 



Finally, the thermodynamic potential r2(/i) — 
£{n{ii)) — fin which is the quantity we shall calculate 
in the grandcanonical ensemble, satisfies 



n{fi) +fi = n{-fi -U)- fi-U 



(7) 



We stress that any approximation method used to solve 
the problem would have to satisfy this symmetry in order 
to yield physically consistent results. 



III. HARTREE + BCS THEORY AT T = 

The starting point of the Hartree -I- BCS theory is the 
Hamiltonian pi). In mean field theory the contribution 
of the interaction term Vk = ~Ucl^ t -k i'^-k.^Ck,-]- to the 
ground state energy can be written as 



(y) = -f/^(c[,^ck.t)(cl 
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The first term is a constant and is the Hartree correction 
to the ground state energy. Note that since the chemical 
potential /i is the derivative of the ground state energy 
w.r.t. n, we can absorb the overall shift of the ground 
state energy due to the Hartree term in /i by adding 
{nU/2) to it. The quantity F^ in the second term is 
self-consistently obtained by minimizing the ground state 
energy w.r.t. i^k and this gives: i^k = Ao/2i?k. The gap 
and number equations in HBCS respectively read. 
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(9) 
(10) 



where i?k = ^^k + ^o- We note that the single particle 
energies ^k in H13CS include a Hartree shift: Ck = fk — 
/i — nU/2, where the last term is the Hartree term. It 
can be easily verified that the HBCS theory satisfies the 
particle-hole constraints on thermodynamics (see Eqs. pi 
[6] and!?]) derived from the attractive Hubbard model. 
Eqs. (pTand ( 10 1 are then self-consistently solved for fi 
and Aq for a given value of n |23 , 24J . The result for the 
chemical potential is plotted in Fig. Q for a given filling. 
The chemical potential is monotonically suppressed as 
a function of coupling. Within the HBCS theory the 
strong coupling expansion of the chemical potential and 
the gap are respectively given by: /iHBCS — ~{U/2) + 



12{tyU){n 

{6tyu)yT 



1) and Ahbcs - {U/2)^l - {1 - n)' - 
(1 — n)^. We next develop a diagrammatic 

formulation to include the effects of quantum fluctuations 

on top of the HBCS theory. 



IV. DIAGRAMMATIC METHOD FOR 

INCLUDING QUANTUM FLUCTUATIONS 

ABOUT HBCS THEORY 

In this section, we outline a diagrammatic approach to 
include Gaussian fluctuations on top of the HBCS the- 
ory. By including quantum fluctuations we expect to 
account for the zero-point motion of the collective mode 
and the virtual scattering of gapped quasiparticles. How- 
ever, since we have already included in our HBCS the- 
ory the leading order Hartree term from the Gaussian 
corrections (see Fig. (II])), we should be careful not to 
double count it. We begin by noting that the Hartree 
term that we have included in the HBCS theory can be 
systematically introduced in the mean field propagators 
using the Luttinger-Ward formalism. The details are out- 
lined in Appendix |X] We next use the Hartree shifted 
propagators to include RPA corrections to the thermo- 
dynamic potential. In order to avoid double counting of 
the Hartree term, we subtract by hand this term from 
the Gaussian thermodynamic potential fig. We have ex- 
plicitly verified that by doing so, we not only restore the 




Hartree 
term 



/\ + +\/ 




n th term 



FIG. 1. RPA loop expansion of the thermodynamic potential 
in the BCS limit. The first diagram is the Hartree term. 



where a is the index for each of the TV flavors. This 
Hamiltonian is invariant under the Sp(27V) symplectic 
group and reduces to the original attractive Hub bard 
model (IT]) after setting A^ = 1. As shown in section VB 



the virtue of working with the above form of interaction, 
where the flavor index a is not conserved, is that it lends 
itself to a systematic expansion in the parameter XjlSi 
around the mean field theory results, exact in the limit 
N -^ oo. Although such an expansion is strictly valid in 
the large-iV limit, it is assumed that the general trends 
of the results found will be correct after setting A^ = 1 
at the end of the calculation. 



correct p-h symmetry for fig , the resulting expression for 
f2g is also rendered manifestly convergent in the absence 
of convergence factors which makes it easier to compute 

However, inspite of the compactness of this approach 
it leads to an unphysical negative compressibility in the 
strong coupling limit (see Fig. 10 ). The result is unphysi- 



A. Particle-Hole constraints for the large-A'^ model 



cal since the system in this limit is comprised of hard-core 
bosons with nearest neighbor repulsion and can hence 
neither collapse (prevented by Pauli exclusion) nor phase- 
separate (ruled out on energetic grounds). The failure of 
the diagrammatic formalism, which was to include quan- 
tum fluctuations on top of HBCS theory, therefore neces- 
sitates a different approach and we next turn to a large- A'^ 
formalism. 



V. LARGE-Af THEORY FOR CROSSOVER ON 
LATTICE 

In this section we give a brief account of the large- A^ 
formalism, which as we shall see, starts with a saddle 
point that is different than HBCS, obeys the p-h con- 
straints appropriate for the large- A^ model at zeroth and 
first order in 1/A^, and most importantly predicts pos- 
itive compressibility for all parameters. In addition to 
satisfying the p-h constraints on the lattice, the way our 
large- A^ theory on the lattice differs from other large- A^ 
approaches in the continuum [TU] is the way we treat the 
fluctuation feedback (see subsection B). Additionally, at 
half-filling there is an emergence of charge density wave 
(CDW) order that the large-A^ theory is unable to cap- 
ture (for reasons discussed later). Hence, we shall work 
away from half-filling where the ground state of the sys- 
tem is a non-degenerate superfluid. 

The starting point of our large- A^ formalism is a gen- 
eralization of the Hamiltonian of Eq. (fTl) to include A^ 
fermion flavors for each spin in the form 



H = -t Y. (4,.Cjo,.+h.c.) 



(11) 



Following the discussion in section ( II B ) , we next de- 



rive a set of p-h constraints appropriate for the large- A" 
model ( 12 ). Using the Hamiltonian ( 12 ) and a p-h trans- 
formation: cl^ ^ = (— l)'rfM,<j, we obtain exact relation- 
ships between thermodynamic variables with n fermions 
(particles) per flavor and ones with 2 — n fermions (holes) 
per flavor. The ground state energy (£"), the chemical 
potential (/i) and the thermodynamic potential (J7) now 
respectively transform as follows: 
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We next develop a functional integral formalism with the 
large- A'^ model and show that it respects the above con- 
straints at zeroth order and also at 0(1/A^). 



B. Functional Integral Formalism 

In this section we shall outline the key steps in for- 
mulating a functional integral approach with the large- A^ 
model. The details are given in AppendixlB] The thermo- 
dynamic properties of the system can be obtained from 
the partition function which can be expressed as a Feyn- 
man path integral over Grassmann fields ^aa and ^qo-. 
We next introduce a Hubbard-Stratonovich field A(x) at 
each X = (xi,r) which couples to I]a ^»at('^)*ia;(T), 
and decouple the quartic fermionic interaction term in 
the action. This makes the functional integral both Gaus- 
sian in the fermionic fields and diagonal in the flavor 
index a. After integrating over these Grassmann vari- 
ables we get an effective action in terms of the Hubbard- 
Stratonovich fields A(a;). It can be easily shown that 
the space- and time-independent saddle point of this ef- 
fective action corresponds to a thermodynamic potential 
that is linear in A'^ and Gaussian fluctuation corrections 
to the saddle point are zeroth order in A^, so that the to- 
tal thermodynamic potential per flavor can be expanded 



as 



n 

N 



= 0n 



>^ 



(15) 



To find the uniform, static saddle point of the effec- 
tive action Sa, we replace A{x) by the space-time inde- 
pendent quantity Aq. The saddle point condition is |27| 
dSo/dAo = 0, which can be rewritten as 



7 4^2£;k 
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(16) 



where -Ek = vCk+^o- "^^^ mean field number equation 
can be obtained from the mean field thermodynamic po- 
tential J7o as 



(dno 



or 
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(17) 



Eqs. (161 and (171 must be solved self-consistently to 



obtain the mean field gap parameter Ag corresponding 
to the mean field chemical potential JL, as well as finding 
the chemical potential which yields the desired density 
n. The results of this calculation are presented as dashed 
lines in Fig. [2] 

It is instructive to show that this mean field theory 
satisfies the particle-hole constraints in the lattice (see 
section n(B)) to the proper order, i.e. to zeroth order 
in 1/A^. From Eq. (pi) we see that this corresponds to 
the chemical potentials on particle and hole sides being 
related by /i(n) — ~'p(2~n). The validity of this equation 
can be seen by replacing /i — >■ — /i without modifying Aq; 
this leaves ( 16 ) unchanged while replacing n — > 2 — n in 



The large U/t limit of this theory can be easily ob- 
tained from the equations. To zeroth order in t/U, the 
chemical potential becomes /I = (1 — n)U/2 and the gap 
parameter is Aq — a/I — (1 — n)^ C//2. We finally em- 
phasize that the essential difference between the previ- 
ously described HBCS theory and this mean field theory 
is the absence of the Hartree term in the latter. Such 
a term, which corresponds to the particle-hole channel 
cannot be easily obtained at the mean field level of any 
functional integral formalism. We recover this important 
contribution in our theory as a 1/A'^ order correction in 
what follows. 

We next expand our action in terms of the fiuctuations 
around the saddle point and truncate to Gaussian order 
or 0{1/N). The Gaussian thermodynamic potential ilg 
can be expressed in terms of the fluctuation propagator 
(see Appendix IbI) , which has poles on the real axis cor- 
responding to the collective mode and branch cuts corre- 
sponding to the two-particle continuum. Using the new 
approximation to the thermodynamic potential per fla- 
vor, Oo(/x) -I- (1/-/V) ^g{n), we next obtain expressions for 
the properties of the system to linear order in 1/A^. At 
this point, we want to emphasize that we do not treat 
the chemical potential fi and the auxiliary field Aq at 



equal footing [55]. Indeed, in our approach the former is 
a thermodynamical variable while the latter is merely a 
parameter in the theory, obtained as the saddle point of a 
variable that is integrated over in the partition function. 
As such, it is not an independent variable but it is defined 
as a function of /i, i.e. Ao(/i) is the saddle point field used 
to calculate the partition function at such a chemical po- 



tential, obtained from the solution of Eq. (BIO I. As we 



make expansions in powers of 1/iV this equation is left 
unchanged, as the saddle point condition is exact to all 
orders 27 . 

In order to calculate the leading order corrections to 
the thermodynamic quantities we next expand the renor- 
malized number equation and the saddle point condition 
using: fj. = fiQ + {l/N)Sfj, and A = Aq -h {1/N)SA. This 
gives us the Gaussian corrections 5/z and SA to the chem- 
ical potential and gap parameter respectively. In order to 
make a connection to the original system with two spin 
components, we set N = 1. 

Finally, we note that, even though our approach to 
the 1/A^ expansion is different from the one introduced 
in reference [TD] it can be shown (see Appendix [e|) that 
the first order corrections to the chemical potential, Sji, 
in both approaches are equivalent while the corrections to 
the gap parameter SAq are different; the latter is due to a 
modification of the gap equation at the 1/N level which 
we do not include. As long as the emphasis is in the 
calculation of the thermodynamics of the system which 
depend solely on the chemical potential, this difference is 
not as relevant. 



Zero-temperature Results for the large- A'^ 
formalism 



Using our formalism we can calculate all thermody- 
namic quantities for the system. In this section we 
present our results, both for the mean field approxima- 
tion and up to linear order in the 1/A'^ expansion; in the 
figures we have set the number of flavors A^ equal to 1 in 
the expansion. 



1. Chemical potential and gap parameter 

The chemical potential and the gap parameter across 
the entire crossover and for a typical density (quarter fill- 
ing, n = 0.5) are plotted in Fig. [21 while the fiuctuations 
are unimportant for small interactions U, the correction 
becomes important at unitarity and in the BEG limit. 
The fluctuations decrease the value of the order parame- 
ter, as well as decrease the value of the chemical poten- 
tial; as we shall see this is related to the Hartree shift in 
the energy of the system. 

We can show that our theory satisfies particle-hole 
symmetry to first order in 1/N. As we can see from 
(14) and the expansion (B7), particle-hole symmetry at 
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FIG. 2. Gap A and chemical potential /x as a function of U jt 
for a filling n = 0.5. The dashed line is the mean field result 
and the solid line is the result which includes fluctuations upto 
order 1/A'^. 



this order implies that 



(18) 



where \i = n{n). This property of ourf2g(/.t) can be 
directly seen from the second line in (B18l; once again 



making the transformation /i — > — // as well as switching 
variables to k — >■ — k and q — > — q we can see the first 
term is left unchanged while in the second one u -^ v. 
Thus, we recover (18). 
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2. Compresstbility 

We next calculate the compressibility of the system 
defined as k = dn/dfj, to order 1/A'^. Using the number 
equation n = —d{no + {l/N)D,g)/dfi, differentiating with 
respect to fi and evaluating the resulting expression at 
fi = JL + 6fi/N, we obtain 



(19) 



Fig. (Is]) shows the compressibility as a function of fill- 
ing in the strong coupling limit {U/t — 30.0) both at 
the mean field level and after including 1/N corrections. 
The positive compressibility in the BEC limit is a check 
that our theory gives physically correct results in this 
regime. The inset in Fig. (|3| shows the behavior of 
compressibility across the crossover. In the strong cou- 
pling limit dn/dfi ex 1/U stemming from the fact that 
/i = (1 — n)U/2 and that fluctuations are unimportant 
for this quantity. 



D. Comparison between HBCS and large-A'^ results 

at r = o 



In this section we shall compare the T = results from 
the large-iV theory with the HBCS results. To begin 




FIG. 3. The compressibility dn/dfi as a function of filling 
n in the strong coupling limit {U/t — 30.0). The dashed 
line is the mean field compressibility. The solid line is the 
compressibility upto 0{1/N) and can be seen to be greater 
than zero and large (but finite) for smaller values of n. The 
inset shows how compressibility changes across the crossover 
for n = 0.5. 
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FIG. 4. Figure shows a comparison of the chemical potentials 
obtained within Hartree shifted BCS theory and large- A'' re- 
spectively. The filling fraction n = 0.5. 



with, we plot the chemical potentials for the two theo- 
ries in Fig. (B. We note that the results from the two 
theories match in the BCS regime; however there is a 
large deviation in the BEC regime. In particular, when 
U/t > 1, 



A'HBCS 



^{U/2) + 12{t'^/U){n-l) 



(20) 



On the other hand, the strong coupling limit of the chem- 
ical potential within the large- A^ theory scales like 



fic^-{U/2){2-n) + 0{tyU), 



(21) 



which has a different leading order term compared to 
the HBCS /i. We already see that there are quantitative 
differences between the HBCS and the large- TV results. 
In order to understand the results (20 1 and (21), we 



next solve the problem exactly in the atomic limit: t/U 



(a) HBCS 
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0.3 

i. 
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(b) large-Af 
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FIG. 5. (a) Compressibility dn/dfi obtained within Hartree 
shifted BCS theory for filling fraction n = 0.5. (b) Compress- 
ibility dn/dn obtained by using large- A'^ theory for filling frac- 
tion n = 0.5. The dashed line is the MF and the solid line is 
the MF-f 1/A'' compressibility. 



and show that the Hartree Mean Field theory gives 
exact answers in this limit. For t/U = the Haniiltonian 
is a single site one 



H = -Un^n^ - fi{n^ + n^) 



(22) 



and has four eigenstates: | 0), jt), \i) and |ti) with re- 
spective energies 0, — ^, — /i and — 2/i — U. To study the 
broken symmetry state we now introduce a fictitious pair- 
ing field h to obtain a Hamiltonian 



H 



^Un^ni - ^i{nf + ni) 



h(c\cl + c.c.) (23) 



After setting up the gap and number, we finally take 
h = 0, solve the gap and number equations for /i and Aq, 
and obtain /x = —U/2 and Aq — U ^n{2 — n)/2. These 
are exactly the values obtained from solving the HBCS 
number and gap equations (|9 10 1 [30 . Physically, the 
result /i = —U/2 in the atomic limit can be explained by 
noting that the chemical potential of the fermions is just 
one-half of the binding energy ~ U for the molecules. 
From these considerations, we conclude that the large- 
N theory gives quantitatively incorrect results for the 
density (ti) dependence of the chemical potential in the 
atomic limit. This also leads to problems with the com- 
pressibility since the compressibility is the derivative of 
/i with respect to n. We turn to this next. 

In order to calculate the compressibility dn/dfi within 
HBCS we note that the Hartree shift to the single particle 
energies can be incorporated into the chemical potential. 
Following |31j we write the renormalized chemical poten- 
tial as fl — fi -\- nU/2, and evaluate dn/dfj, as 



dn dn/dfl 

d4i ~ 1 - {U/2)dn/djl 



(24) 



The quantity dn/dfi can be calculated from Eqs. ([9] 
and PI 



dn 



A, 



2 V^ ^ 
k ^k 






(25) 



Fig. (5) shows a comparison of dn/dfj, as obtained from 
HBCS using Eqs. (24 25) and large- iV respectively. In 



the strong coupling limit the HBCS compressibility scales 
as 



dn 
d^i 



U 
12t2 



1 



24t2 



(26) 



This is again expected on general grounds, when one 
notes that the chemical potential /i, written in powers 
of t/U, has a zeroth term equal to —U/2 (atomic limit). 
Any n dependence of fi is hence at least 0{t^/U), which 
implies that the compressibility must increase with U. 
From these considerations we come to the conclusion that 
the compressibility should be a monotonically increasing 
function of U/t. In contrast, the compressibility for large- 
N scales like dn/dji ~ 2/U. 

To summarize we find that the much simpler BCS plus 
Hartree theory works better in the BEC limit compared 
with the more sophisticated large- A^ approach where we 
included the 1/A^ Gaussian fluctuation corrections to the 
saddle point result, with both approximations satisfy- 
ing the particle-hole constraints on the thermodynamics. 
By better we mean that BCS -I- Hartree reproduces the 
atomic limit behavior of the chemical potential and the 
strong coupling behavior of the compressibility expected 
for a BEC of hard-core lattice bosons, while the large- 
A^ approach (with A^ set equal to 1 at the end) does 
not. We can also compare our results with the available 
Quantum Monte Carlo data, which however is only for 
the two-dimensional attractive Hubbard model. We find 
that the results for the chemical potential [32] and for the 
compressibility [3T1 at moderately large \U \ /t are both 
in good (semi-quantitative) agreement with the Hartree 
+ BCS theory. Although there exists no QMC data on 
the 3D attractive Hubbard model at T = 0, we note that 
fiuctuations should be less important in 3D than in 2D. 
It is therefore reasonable to expect that the agreement 
between Hartree -I- BCS theory and QMC should only 
improve in 3D. 

The comparison between the two approaches (large- 
A^ and HBCS) is quite surprising and unexpected. We 
should emphasize that the two theories start with quite 
different mean field solutions (or saddle points). The 
BCS -I- Hartree solution incorporates both the particle- 
particle (p-p), or pairing, and the particle-hole (p-h) 
Hartree physics on an equal footing at the mean field 
level. The large-A^ solution, on the other hand, is de- 
signed to focus only on the p-p channel at the saddle 
point level, and include all other effects as fluctuations 
about the saddle point. One might have thought that 
since the Hartree correction to the thermodynamics is in- 
cluded at the 1/A^ level (along with higher order terms), 
the large-A^ approach would "go beyond" the simpler 
BCS -|- Hartree approach. But we find that "more" (dia- 
grams, for instance) is not necessarily "better" for quan- 
tum many-body systems! 

It might also be worth contrasting the optical lattice 
calculations presented here from BCS-BEC crossover in 
the continuum. In the continuum, one does not in gen- 
eral have a Hartree term in thermodynamics, which is 



proportional to both the interaction and the density (ex- 
cept in the BCS hmit). The reason is as follows: the bare 
interaction g{A) actually goes to zero as the ultraviolet 
cutoff (inverse range of potential) goes to infinity. Thus 
a "bare Hartree term" proportional to g{A)n vanishes 
throughout the crossover. Also there can be no term 
proportional to the renornialized interaction a^ in the 
ground state energy in general, since that would diverge 
at unitarity! As shown in Ref. [S], the Hartree diagram 
in the Gaussian fluctuation correction to the BCS theory 
does indeed lead to the expected Hartree correction of 
relative order kfUg in the BCS limit. But it is not pos- 
sible to isolate "the Hartree correction" to the ground 
state energy or chemical potential at arbitrary values of 
l/kfttg in the continuum problem. 

In the following section we compute the superfluid den- 
sity on a lattice at T = and show that because of the 
broken translational invariance, the superfluid density is 
not equal to the total density as is the case with the 
continuum. 



VI. SUPERFLUID DENSITY 

It has been shown in Ref. [23, on quite general 
grounds, that the superfluid density of a translationally 
invariant superfluid possessing time reversal invariance at 
T = 0, is equal to the total density. For a one component 
system, barring pathologies (e.g. He^-He^ mixtures), the 
statement can be proved using the Gibbs-Duhem relation 
and Landau's two-fluid model. A phase twist put in the 
boundary conditions for the order parameter in a transla- 
tionally invariant system is uniformly distributed across 
the system. The situation is different on a real optical 
lattice - because of the broken translation symmetry, the 
many-body wavefunction has a very small amplitude be- 
tween the lattice sites and it is energetically advantageous 
to distribute the phase twists at these locations. As a re- 
sult the superfluid density, which is the response of the 
system to this phase twist, turns out to be different on a 
lattice. Indeed one can show, using Kubo formalism, that 
for a translationally invariant system the paramagnetic 
part of the current-current correlation function vanishes 
due to the commutativity of the total momentum op- 
erator with the Hamiltonian. The superfluid density in 
such a system is therefore entirely given by the diamag- 
netic part of the response and turns out to be equal to 
the total density. In a discrete lattice model, the total 
momentum operator does not commute with the Hamil- 



tonian and hence the paramagnetic part of the response 
is non-zero. Consequently, the superfluid density differs 
from the total density on a lattice even at T = 0. 

The superfluid density is computed by comparing the 
free energy F{n) = il + iin oi the gas at rest with the 
free energy of a gas moving with a superfluid velocity 
Vs — Q/(2m) in the limit Q — J' 0; indeed, F{Q,n) — 
F(0,n) = ^Usmv^^ so that ^5] 



/d^F 



(27) 



We can relate this derivative of F with derivatives of il 
recognizing that F{Q, n) = ^{Q, l^{Q, n))+ii{Q, n)n and 
thus 



d^F 



dQ-' 




where we have used the number equation at Q = in the 
last line. 

Following [53] we calculate, within our large-A^ for- 
malism, the superfluid density on a lattice as a response 
of the system to a phase twist on the order parameter. 
We give the details of this calculation in Appendix [F| and 
present the results here. 

The mean field superfluid density for the large- TV the- 
ory is given by 



1 

E 



d0^ 



1.1, A 



1- 



Ek 



cos{kxa) 



(29) 
(30) 



From this formula we can see that in general we have 
n^ < n. The Galilean invariant result in the dilute, con- 
tinuum limit is recovered in the limit of small density 
ji and interaction U/t, in which the chemical potential 
/i is near the bottom of the band. Thus, for the mo- 
menta contributing to the sum we have cos(fc2;a) ~ 1 and 
n^ — n. 

Next, we obtain the 1/N corrections to ps by includ- 
ing Gaussian fluctuations in the calculation of thermody- 
namic potential (see Appendix If] for details). We obtain 
the following expansion for the superfluid density up to 
Oil/N): 



Us ^ n,+ 



N 



d fd^no 



d^i \ 86'^ 



6fi 






mAo 



dSg/P 



-7^ ) (2«(m)) 



(31) 



where d/dp. = d/dfi+{dn/dAo)d/dAo and a(/x) is given by 



v{fi) = - 



3 Euid'^u/dkl)C^/El 
8Aoa2 E^il/El) 



(32) 
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FIG. 6. The superfluid number density Us as a function of 
U /t for n — 0.5. Tiie dasiied line shows the mean field su- 
perfluid number density, the solid line is the result when 1/N 
corrections are included, and the dotted line is l.Olt/U show- 
ing that in the strong coupling regime the superfluid density 
scales like t^ /U . 



We finally set iV = 1 and plot the mean field super- 
fluid density and the one including 1/N corrections as a 
function of coupling strength in Fig. ^. As it can be 
seen, Us falls off like t^ /U in the strong coupling limit 
which can be explained by noting that the system in 
this regime comprises of hard-core bosons on a lattice 
with an effective effective mass ~ U/t^. In other words, 
the tightly bound pairs in the BEC limit can hop only 
through virtual ionization and hence has an hopping pa- 
rameter ^ t^ /U. Further, Gaussian fluctuations reduce 
the superfluid number density from its mean field value 
across the entire crossover with an increased suppression 
in the strong coupling regime. This is expected because 
the BCS mean field theory reduces the problem to one of 
non-interacting Bogoliubov quasiparticles with a gapped 
excitation spectrum. However, the low lying excitations 
in the strong coupling regime, are the gapless collective 
modes of the composite bosons, which are not captured 
by the BCS mean field theory. 

In the next section we calculate the critical tempera- 
ture within the large- TV theory and show that the large- 
N theory for the attractive Hubbard model, inspite of its 
above limitations, predicts the correct qualitative trends 
for Tr in the two limits. 



VII. CRITICAL TEMPERATURE 

Let us finally calculate the critical temperature Tc of 
the fermionic gas in the lattice, as well as the pairing 
temperature T* . Just like in the continuum fF, these 
two temperatures are approximately the same in the BCS 
limit and widely diverging in the BEC limit. In the for- 
mer, the formation of Cooper pairs and their condensa- 
tion are governed by the same physics. In the BEC limit, 
on the other hand, the temperature for the formation of 
pairs will be of the of the order of the binding energy 



(which is proportional to L'^) while the critical temper- 
ature will be decreasing as U increases, as the effective 
mass of the pairs increases, as we shall show. 

To calculate both T* and Tc we need to consider the 
temperature at which the t-matrix has a divergence at 
zero energy and momentum for a given chemical potential 
/J,. This condition for the inverse temperature /? = T~^ 
is 



i.5:^tanh( 



/^a 



U 



) 



(33) 



The two temperatures differ, however, in the equation 
of state that is used to calculate the density. The pair- 
ing temperature is obtain using the mean field approxi- 
mation to the thermodynamical potential, in which only 
fermionic excitations are included (i.e. pair breaking at 
that temperature). The calculation of the critical tem- 
perature corresponds to the addition of the effects of 
Gaussian fluctuations in which bosonic excitations (Gold- 
stone modes) are included, leading to a large renormal- 
ization of the equation of state. In order to set up the 
number equation we use the same functional integral for- 
malism that we developed at zero temperature. The 
mean field thermodynamic potential at high tempera- 
tures /3 < /3c is: 



^0 



/? 



Y^Hi 



-/3«k 



(34) 



The mean field number equation at T — Tc is then given 
by dflo/d^ = — n, or 



E 



exp(^a) + 1 ■ 



(35) 



(33 35) self-consistently we obtain the 
for a given density n and denote it by 



Solving Eqs. 
mean field /3c 

po = (r*)-i. 

The calculation of the Gaussian contribution to the 
thermodynamic potential at these high temperatures (at 
which the gas is normal) is also similar to the one pre- 
sented at T = setting A = everywhere. Thus, the 
fluctuation propagator M is diagonal, with 'M.ii{q) = 

r(g) = l/U + Ek(l - / - /')/(*« -i-O where 
iqi = i2TTl//3 are Bose Matsubara frequencies and /(^) = 
l/[exp(/3^) + l] is the Fermi distribution function. Hence, 



f^, = ^5]lnr(g) 



(36) 



Using our large- A^ formalism, we shall obtain the 1/N 
expansion of the diagrammatic approach to this prob- 
lem, which was discussed by Nozieres and Schmitt-Rink 
(NSR) 3J. Following NSR, we maintain the same form 
of the t-matrix equation (33). Nevertheless, just like at 



T = the 1/A^ corrections to the thermodynamic po- 
tential will renormalize the chemical potential and hence 
change /3c. The inverse temperature (3 and the chemical 
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FIG. 7. The transition temperature Tc/t as a function of U /t 
for n = 0.5. The dashed line is the mean field transition tem- 
perature denoted by T* /t. The solid line includes fluctuations 
upto 0{l/N). 



potential /i are in general independent thermodynamic 
variables but are related to each other aXT — Tc through 
the t-matrix equation. Expanding Eq. (33) upto 0(1/Af) 
we obtain 



dn^ 1 f d'^Vto 



9Ao N V^Aoa// 



5^ 



dAodp 



6(3 







(37) 



In the above equation all the derivatives are evaluated at 
A — 0, ij, — JI and (3-/3^. Setting the coefficient of the 
1/N term to zero we obtain 



6(3^ 



{d^n^/dA^dtj) 



Sjj, 



id^no/dAod/3) 
Similarly we expand the number equation and obtain 



(38) 



ana 

9/i 



1 

N 



ffHh 



(5/i 






513- 



dflg 



(39) 



and setting the coefficient of the 1/iV term to zero we get 
(dng/dfi) + (d^no/dfidp) (5/3 



(5/.J = — - 



(d^Qo/dtx^) 



(40) 



Eqs. (38 40) are then simultaneously solved to obtain 
the 1/iV corrections to /?° and JI. We obtain the mean 
field temperature scale T* and the critical temperature 
Tc after including 1/A^ corrections from T* = 1//3J? and 
Tc = il3c + <5/3)~^ (after setting A^ = 1) respectively. 

The results are shown in Fig. (It]). As expected, there 
is a large deviation of Tc from its mean field value T* 
in the strong coupling regime. The phase diagram is as 
follows: above T* the system is a normal fermi gas, for 
temperatures below T* and above Tc there are preformed 
uncondensed pairs and below Tc we have a condensate of 
pairs. In the weak coupling limit Tc approaches the BCS 
value Ao/1.75. Beyond the BCS regime, the pairing tem- 
perature grows linearly with U/t while Tc goes through 
a maximum near unitarity and then falls off as t^ /U . 



FIG. 8. This plot shows the scaling of transition temperature 
Tc/t with coupling U/t {U is on-site interaction and t is hop- 
ping matrix element) for a filling of n = 0.5 in the two limits. 
The solid black line is Tc/t, the dotted line is the zero tem- 
perature gap parameter Ao/t (rescaled by a factor a = 0.57), 
the dashed line is zero temperature superfluid stiffness Ds 
(rescaled by a factor 7 — 6.67). One can see that Tc scales 
like Ao in the BCS limit and like n^ in the BEC limit. 



The scalings of Tc in the two regimes are summarized 
in Fig. ([8]). In the weak coupling regime the pair break- 
ing energy scale Aq is much smaller than the energy scale 
set by the superfluid stiffness, D^ = Ust and hence the 
transition temperature Tc is governed by the zero tem- 
perature gap. In the strong coupling regime, the energy 
scale for phase fluctuations is the smaller one compared 
to the pair breaking one and hence the scale for Tc is dic- 
tated by the zero temperature superfluid stiffness. The 
precise value of U/t for which Tc goes to a maximum 
depends on filling (see Fig. l9|. As explained earlier, 
the reason why Tc scales like i /C/ in the BEC limit can 
be explained by considering hardcore bosons on a lat- 
tice. A simple second order perturbation theory in t/U 
then gives an effective hopping parameter proportional 
to t^ /U for the composite bosons. Since the superfluid 
stiffness is proportional to the effective hopping param- 
eter and Tc in this regime is governed by phase fluctua- 
tions of the lattice Bose gas, therefore Tc ^ t^ /U in this 
limit [161 I34j . Lastly, we have compared our large- A^ re- 
sults with a Hartree shifted NSR theory (HNSR), which 
uses Hartree shifted propagators in the scattering matrix. 
The details are given in Appendix [G] We find reasonable 
agreement between the two approaches at BCS and BEC 
limits (see Fig. 



12 1. However, there are quantitative dif- 



ferences between the two results around unitarity which 
we do not understand at this stage. 



VIII. CONCLUSIONS 

In this paper we have addressed the BCS-BEC 
crossover on a 3D optical lattice. We have developed a 
simple Hartree -I- BCS theory that satisfies the p-h con- 
straints imposed by the lattice on the thermodynamics. 
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FIG. 9. (color online) Plot shows Tc/t as a function of U /t 
for various fillings: n = 0.4 (black), n = 0.5 (red), n = 0.6 
(green), n = 0.7 (blue) 



due to pair breaking and in the strong coupling the su- 
perfluid order is destroyed by phase fluctuations of the 
lattice Bose gas. 
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Appendix A: Diagrammatic Approach for the 
attractive Hubbard model 



Since, inclusion of Gaussian fluctuations on top of the 
HBCS theory led us to an unphysical negative compress- 
ibility in the strong coupling limit, we were forced to start 
from a different saddle point. We developed a large- A^ 
approach for the attractive Hubbard model in 3D, where 
the large- iV saddle point did not include the Hartree shift 
but still respected the respective p-h constraints of the 
large- iV model. Most importantly, inclusion of Gaussian 
fluctuations led to a finite and positive compressibility 
for all parameters. We calculated the ground state chem- 
ical potential, gap, ground state energy etc. away from 
half-filling. The superfluid density at T = on the lat- 
tice was found to deviate from the total density and in 
the BEC limit was determined by the single-boson hop- 
ping matrix element which scales as t/U . However, we 
flnd that the large-A theory predicts quantitatively in- 
accurate results for the chemical potential in the strong 
coupling limit, and qualitatively incorrect trends for the 
compressibility across the crossover. A comparison with 
the HBCS theory, which correctly recovers the atomic 
limit and predicts the right qualitative trends for com- 
pressibility, reveals that the large-A theory on the lat- 
tice, although considers a larger number of diagrams, is 
in fact inferior to the simpler Hartree shifted BCS the- 
ory. The limitation of the large- A" approach is explained 
by noting (i) the importance of Hartree shift in lattice 
problems, and (ii) inability of the large-A^ approach to 
treat particle-particle and particle-hole channels at equal 
footing at the saddle point level. 

Inspite of the limitations of the large- A^ approach in 
describing the two-component fermionic system on the 
lattice at T = 0, we obtain correct trends for the critical 
temperature within this theory by approaching the super- 
fluid state from above T^. The transition temperature is 
shown to scale like two different ground state quantities 
in the two regimes: in the BCS regime T^ scales like the 
gap, while in the BEC regime T^ scales like the zero tem- 
perature superfluid stiffness. These two different scalings 
show that in the weak coupling regime coherence is lost 



In this appendix we develop a diagrammatic formula- 
tion of the crossover problem in the lattice in order to in- 
clude gaussian fluctuations on top of the Hartree + BCS 
theory. The starting point of this discussion is the attrac- 
tive Hubbard Hamiltonian (fTj) . The thermodynamics un- 
der this Hamiltonian obeys the p-h constraints discussed 
in Eqs. ([s] |6J and[7|). It is easy to see that starting with 
Hamiltonian ([l| if we were to develop a functional inte- 
gral formalism (like the one used for the large A^ model) , 
we would flnd that the saddle point violates the p-h con- 
straints. The reason can be traced back to the choice of 
Hubbard-Stratonovich fleld (either p-h or p-p channel). 
At this stage let us anticipate that a Hartree shift to the 
mean fleld chemical potential would correct this problem. 
We use a Luttinger and Ward formalism [35J to system- 
atically introduce a Hartree shift in our single particle 
propagators. Introducing the Luttinger- Ward functional 
$[0] we write the thermodynamic potential il as 



f2 = $ + Tr In G - Tr EG 



(Al) 



where Tr = (1//3) ^j, ^^ tr and the self-energy S is ob- 
tained by evaluating the functional derivate of ^[G] at 
the exact Green's function 



.^..,o,.. 



(A2) 



Note that the relation 5][G] = S is independent from 



the Dyson equation G 



S, where Go is the 



non-interacting Green's function. In the Luttinger- Ward 
formalism $ is obtained by summing up an inflnite se- 
ries of closed diagrams without any self-energy insertion 
(generally called skeleton diagrams) and replacing all free 
propagators by fully interacting ones. At the mean field 
level we need to retain only the first diagram in the series 
and thus 

$[G] = -f/(TrGi2)(TrG2i) + [/(Tr Gn)(TrG22) (A3) 

We define (5$[G]/(5Gii = UTTG22 = S which implies 
S<^[G]/SG22 = UTrGii = -S. We can further asso- 
ciate TrG2i — TrGi2 with the Hubbard-Stratanovich 
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field A and therefore S<^[G]/SGi2 = -[/TrGzi = -A. 
The Luttinger-Ward functional at the mean field level is 
therefore given by 



$fGl = - 



A2 



U U 
and the self-energy matrix is given by 



S -A 
-A -S 



(A4) 



(A5) 



We next use the form of the free Green's function given 

by 



'o' - [ " 



Q-l ^ I 'l-l^n - Ck + A^ 










(A6) 



and the Dyson equation to calculate the inverse of the 
full Green's function 






Ck + M 

A 



Ek - M 



(A7) 



Note, the Green's function in Eq. (A7| has its single 
particle propagators Hartree shifted. 



1. Mean field theory at T = 



Using Eq. (Al) and the fact that TrS = we can 



obtain an expression for the mean field thermodynamic 
potential 



no 



u 



'r-Y.iE.-i. 



E2 
U 



(A8) 



where -Ek = ^/i"^ + Aq and ^k = ek — A* + S. This form 
of thermodynamic potential, as anticipated earlier, obeys 
the correct particle-hole constraints. Then the spatially 
uniform, static saddle point at T = is given by the 
following condition 



SAo 



or 






U 



(A9) 



The mean field number equation can be obtained from 
the condition 



dVLa 



-n or n 



El 



and the Hartree shift E is given by 



dT, 



or E 



El 






in 

Ek 



(AlO) 



(All) 



Eqs. (A9), (All), and (AlO) arc then solved self- 



consistently and we obtain the mean field values for Ag, 
fi and E. 



2. Gaussian fluctuations at T = 

In order to go beyond the mean field approximation 
we next consider fluctuations of the order parameter A 
around its static saddle point value Ag and expand the 
action Sa to Gaussian order. The first order term van- 



ishes due to the saddle point condition ( A9 ) and we ob- 
tain 



Sa — Sq + Sn 



(A12) 



The mean field piece 5*0 has been defined above and 
Gaussian piece has the form 



Sg^lj2iv*iq)v{-q))M{q)f ^^'^^ 



v*i-q) 



(A13) 



where iqi — i2Til/ j3 are the Bose-Matsubara frequencies 
and the matrix elements of the inverse fluctuation prop- 
agator M are given by 



Mn(g)=M22(-(7)-l 



ik,-,,M. 



G",,ik)GUk + q) 



1 

U 



E 

k 



4"k' 



<^lvl 



iqi - -Bk - -Bk' iqi + Sk + -Bk' 



(A14) 



and 



Mi2((j) = M2i(<?) = Y. G?2(fc)G?2(fc + q) 



2/c^,k 



E 



UkWk'WkWk' 



1 



iqi + Bk + Bk' iqi - E] 



'k — -C-k' / 



Here G is the same Nambu propagator defined in Eq. 



(|A7|) with A = Ao,ul = l-vl = (1/2)(1 + ^k/B^k) 
are the standard BCS coherence factors and k' = k + q. 
While calculating the thermodynamic potential includ- 
ing Gaussian fluctuations we need to remember that the 
first term in the Gaussian part (J7g) is indeed the Hartree 
term (— E^/t/). Since the Hartree contribution has al- 
ready been included at the mean field level to preserve 
particle-hole symmetry, we need to take it out from il^ 
to avoid double counting. Writing the partition function 
upto Gaussian order and integrating out the Gaussian 
fiuctuations we obtain the Gaussian contribution to the 
thermodynamic potential 



"" = ^El^^f^DetM(,))e 



2/3 



g^.q 



22 



tq„a+ 



U 



(A16) 



where the matrix elements Mn etc. have been rescaled 
as Mil ~> C/Mii. It is easy to see that in the limit of 
large iqi 



DetM(q) - 1- 



t/2 



{^qly 



Y^/ ,2„,2 



. k 



^k"k' 



2 2 
"k^'k' 



(A17) 
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However 



In 



M22 



In 



2U 
iqi 



E(^ 



MuMu/ — V\,V 



'k"k 



ki^k' 



(A18) 



and hence the Matsubara sum ^-^ hiMn Det M/M22 
without the convergence factor diverges for large iqi. 
However, the correct Jig also has a correction term given 



by EVt/ 
U 



2«E^E [Gii(k)G22(k') + Gn(k')G22(k)] 



^^[(M22-l)e-^'°^+(Mn-l) 



ViqiQ^ 



2p 



(A19) 



Upon changing the sign of q in the second term of the 
second line and noting that the sum over q is over both 
positive and negative values we have for large iqi 



U 



^E??E("k"^'-«^kOe+-°" (A20) 



213 ^iqi 

l(Ji,q k 



which exactly cancels the linear term in the large (iqi) 
expansion in Eq. (AI8). 



Summing up the above results we obtain the Gaussian 
correction to the thermodynamic potential 



n„ 



1 

2p 



El- 



«9»i,q 



(Mu(q)/M22(g))DetM(g) 



exp [t^Ek {«'l{ili - ^k - -BkO - vlvlliiqi + Sk + E^,)}] 



(A21) 



where we are justified to drop the converg ence factor 
g+MjiO fj-om the right hand side of Eq. (A21) since 



in the large {iqi) limit, the leading order term in the 
sum is now of the order l/{iqi)^ and thus the Matsubara 
sum is convergent. Thus the same scheme that restores 
the correct particle-hole symmetry in our theory, also 
makes the Matsubara sum convergent at the Gaussian 
level. To evaluate the Matsubara sum in Eq. (A21) we 



analytically continue in the complex plane and convert 
the sum over the bosonic Matsubara frequencies to an 
integral over a closed contour enclosing the imaginary 
axis counter clockwise ^^ — > §{dz/2'Ki)nB{z) where 
nsiz) is the Bose distribution function. We evaluate 
the integral over z along a contour parallel to the Mat- 
subara axis: z —> 0^ + iy keeping in mind that the 
phase of lnMii(q, j/)/M22(q, J/) and the imaginary part 
of (Mii(q, y) — 1) are both odd functions of y and hence 
do not contribute when integrated over positive and neg- 
ative values of y. Therefore, we obtain at T = 



n„ 



«<?!,q 



Mil Det M 



— Vln, 

2/3 f-^ VM22exp(-S7C/) 



/•OO 

= / dy/(27r) V[ln(DetM(2/))-2Re(Mii-liJ|\22) 

To obtain Aq , /i and E including gaussian corrections we 
start with a grand canonical ensemble and treat both /i 
and E as thermodynamic variables. For convinience we 
switch to /i = /i — E and E as our independent variables. 
Then, the thermodynamic potential can be written as 
ri(/2, A(/i); E) = A{fl, A{p,)) + E^/f/, where the function 



A(il,A(fl)) has no explicit dependence on E. The gap 
A(/i) is obtained from the saddle point Eq. (A9|. To 
obtain the number equation and the equation for E we 
construct a function F{jj., E) = ^{p., E) + (/i + E)n. The 
condition for E is then given by 



rOF 
V5S 



= or E = -- 



nU 



(A23) 



The number equation reads 



dF 



or 



dA 



dA 
dA 



dA 
djl 



(A24) 



We next switch to a canonical ensemble and for a fixed 
n numerically calculate A[jl,A{jl)\ = Aq[^,A{^)\ + 
A,[/i, A(/i)]. Eq. (A24) then gives the value of the renor- 



malized Hartree shifted chemical potential fi for the cor- 
responding value of n which when combined with Eq. 
( A23 1 gives the renormalized chemical potential // with- 
out the Hartree shift. The problem with this diagram- 
matic approach is that it predicts an unphysical negative 
compressibility in the BEG limit. In Fig. (fTol) we have 
plotted /x as a function of n for U /t — 20.0. Glearly, the 
slope of /I versus n is negative for a large range of n indi- 
cating negative compressibility. However, we know that 
in this limit the system is a lattice Bose-gas with a hard- 
core repulsion coming from Pauli exclusion and a nearest 
neighbor repulsion proportional to t^ /U. Hence the sys- 
tem is stable in the BEG limit and the negative compress- 
ibility within the diagrammatic approach is therefore an 
unphysical result. 
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i-14 




-16^ 



FIG. 10. The chemical potential /i plotted as a function of 
the filling n for U/t — 20.0. Note that the slope is negative 
upto n ~ 0.7 indicating a negative compressibility. The range 
of fillings for which dn/d^ < increases with U/t, so that 
eventually for very large couplings the system is unstable for 
all fillings. 



Appendix B: Details of large- A'^ formalism 

The thermodynaniic properties of the system can 
be obtained from the partition function in the grand- 
canonical ensemble Z(^,/3), where /3~^ is the tempera- 
ture T of the system. Indeed, Z is related to the ther- 
modynamic potential as il{^, /3) = — f3^^ In Z. This par- 
tition function can be expressed as a Feynman path in- 
tegral over Grassmann fields ^aa and "^aa 



Z= D^^„D^^,eM-S^) (Bl) 

with the action in imaginary time r 

(B2) 
The quartic fermionic interaction term in the Hamil- 
tonian can be decoupled by introducing a Hubbard- 
Stratonovich field A (a;) at each x = (x^, r) which couples 
to ^^ ^ ia^i'r)'^ ia^ij) . The partition function can then 
be written as Z = J DADA*D^^aaD'i'xacrexp{-S^Sf,A) 
with a full action 



5*, A = Mr y^ ( 



N I A{x) 
U 



(B3) 



function G- {t,t') is given by 

{-dr + fJ.)d^,j + tS^i,jy A{x)S^,j 

A*{x)Sij {-dr - fi)Sij - i(5<ij> 

X(5(t-t'04) 

with the notation S<:ij^ — 1 only if the i and j sites 
are nearest neighbors and zero otherwise. The functional 
integral is now both Gaussian in the fermionic fields and 
diagonal in the fiavor index a. After integrating over 
these Grassmann variables we get 



Z= f DADA* cxp{-Sa) 



(B5) 



with an effective action Sa which only depends on the 
auxiliary fields A(a;) in the form 



Sa^N dx 



A(x) 
U 



'N 



f dxTT\nG-^[A{x)] (B6) 



where J dx = J2i I ^t. 

Assuming that the saddle-point auxiliary field is space- 
and time-independent (i.e. A(x) — Ao), the thermody- 
namical potential Vt is of the form ^{n, /3) ~ iVfJo = 
Sa{A{x) — Ao)//3. Fluctuations around the saddle point 
yield corrections that are smaller than this term by pow- 
ers of 1/A'^; thus the full thermodynamic potential will 
be expanded in the form 



N 



= 0n 



>^ 



(B7) 



1. Saddle point approximation - Mean field theory 
at T=0 

To find the uniform, static saddle point of the effective 
action Sa, we replace A(a:) by the space-time indepen- 
dent quantity Aq. Fourier transforming all the fields to 
the reciprocal (momentum) lattice and Matsubara fre- 
quencies, the effective action is given by 



^a[Ao]=A^ 



U 



N 



E 



trlnG(7^(/c)EEiVS'o (B8) 



with 



Go'(fc) 



*fcn ^ Ck Aq 

Ao ikn + Ck 



(B9) 



where zfc„ = {2n + l)TTi/(3 are the fermionic Matsub- 
ara frequencies. The saddle point condition is [27] 
dSo/dAo — 0, which can be rewritten as 



1 

U 



y— 



(BIO) 



' fdr'Yl V^L(^)G-^(r, r')V,,a(r')) , where i?k = v^ 

•^ ija in the mean-field 



Aq. The thermodynamic potential 



where we have introduced the Nambu spinors tpjai''') — 
(^iQt ('''): ^ia^ (■''))■ The inverse Nambu-Gorkov Green's 



in the mean-field approximation is then 

^0 = So/P = ^ - y{E^ - Ck) 

k 



(Bll) 
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The mean field number equation can be obtained from 



<9m Jt,v V V ^k 



(B12) 



Eqs. (BIO B12) must be solved self-consistently to ob- 



tain the mean field gap parameter Aq corresponding to 
the mean field chemical potential /I, as well as finding the 
chemical potential which yields the desired density n. 



2. Gaussian fluctuations at T = 

In order to go beyond the mean field approximation we 
must consider perturbations of the auxiliary field A(a;) 
beyond the saddle-point, in the form 



A(a;) = Ao + r/(a;) 



(B13) 



where the complex bosonic field ri{x) describes space- 
time dependent fluctuations around the uniform static 



value Aq. We next expand the action Sa in Eq- (B6) 



to quadratic order in rj, using that the saddle point con- 
dition (BIO I ensures that there is no term linear in rj. 
Thus, the action is of the form Sa — NSq + Sg + ... with 
a Gaussian piece of the form 



Sa^l^J2('^*ilW-q))Miq) 



m ,q 



v*i-q) 



(B14) 



where iqi = i2TTl/(3 are the Bose-Matsubara frequencies 
and the matrix elements of the inverse fluctuation prop- 
agator M are given by 

Mn(g) = M22(-g) = 77 + E G^2(fc)G?i(fc + qjElb) 

ikn ,k 



u 



1 

u 



E 



2, ,2 



utu 



k^k' 



2„,2 



vtv 



k'^k' 



m 



Elr — E\, 



iqi + £^k + E^' 



and 

Mi2(9) = M2i(9) = Y. G?2(fc)G?2(fc + q) 

ik„,k 
^ , 1 1 

k 



(B16) 



m 



Ev + -Ek 



iqi 



Eu 



Here we use the standard BCS notation w^ = 1 — ti^ = 
(l/2)(l+a/Ek)andk' = k + q. 
Writing the partition function upto Gaussian order 

Z ~ exp(-A^S'o) / DrjDr]^ exp(-S'g) (B17) 

and integrating out the Gaussian fluctuations we obtain 
(see Appendix [C] for details) the Gaussian contribution 
to the thermodynamic potential 



^. = ^Eln(^DetM(,)]e-- 



•■22 



(B18) 






1 



ln(C/2DetM(q)) 



-E 

2 ^ 

k 



(ul-vl). 




10 15 

q(l,l,l) 

FIG. 11. Spectrum of excitations which contribute to the 
leading 1/N corrections to the thermodynamic potential plot- 
ted at unitarity for n = 0.5 and along the main diagonal 
{q-x — Qy = Qz) of the Brillouin zone of a 50 x 50 x 50 lat- 
tice. The solid line is the collective sound mode given by 
Det M = 0. The shaded region is the two-particle continuum 
given by the branch cut of the fluctuation propagator. 



In a previous article [6] some of us showed that this 
Gaussian fluctuation contribution can be physically in- 
terpreted by analytically continuing the bosonic Mat- 
subara frequency to the real axis iqi ^ z = ui + zO+. 
We are thus led to the study of the analytical properties 
of lnDetM(q, z). The zeroes of DetM(q, z = a;o(q)) 
(which correspond to poles of the fluctuation propagator 
M^^) correspond to the frequencies wo(q) of collective 
excitations of the system with momentum q. These ex- 
citations are the q — ?► Goldstone modes of the order 
parameter in the broken symmetry superfluid state. Ad- 
ditionally, the fluctuation propagator has branch cuts on 
the real axis originating at Ec{q) = ±min(i?k + £^k+q)- 
These branch cuts represent the two-particle continuum 
of states for scattering of gapped quasiparticles. The 
Gaussian contribution (B18) can be then rewritten as 



0„ = - 
^ 2 



E 



wo(q) - E,,{q) - 



-^^('i) diU 



(5(q,a;) +71 



(B19) 
where the last integral describes the contribution of the 
virtual scattering of quasiparticles with a phase shift 
(5(q, w) whose particle continuum begins at Edq) and 
the last term TZ comes from using the correct conver- 
gence factors in the calculation (see Appendix B). 

To illustrate this excitation spectrum we plot in Fig. [TT] 
the two particle continuum and the collective excitations 
along the main diagonal q (1, 1, 1) of the Brillouin zone, 
at unitarity and for n = 0.5. For small q, the collective 
excitation spectrum is linear indicative of sound modes, 
eventually hitting the two-particle continuum. 

In the BCS limit, the contribution of the collective 
mode is negligible due to phase space restrictions and the 
two-particle continuum dominates. In the BEG limit, the 
two-particle continuum lies at a much higher energy scale 
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and the low-energy excitations are entirely given by the 
gapless sound modes. Further, at half-filling, one would 
expect the collective excitation spectrum to be gapless 
at q = (tTjTTjTt) indicating new Goldstone modes due to 
the onset of CDW order [25j. However, since we only 
decouple the quartic interaction in the p-p channel, we 
do not see the CDW order and hence there is no softening 
of (tt, tt, tt) mode at half-filling within our theory. 



3. Corrections of order 1/N 

In order to calculate the leading order corrections to 
the thermodynamical quantities, such as the chemical po- 
tential in this case, we write it as the expansion 

Naturally, given that the gap parameter Aq is a function 
of /^, it will also have an expansion in powers of 1/A^ 
derived from this expansion. Indeed, 

Ao = A;;+^(5Ao; (B21) 



M = Ai 



(B20) 



<5Ao 



dAo 



N 



6fi 



(d^So/dfidAo) 



Next, expanding the number equation to linear order 
in 1/A'^ and remembering that in calculating derivatives 
with respect to ^ (which we denote here as d/d/i) the 
parameter Aq actually changes with n, we get 



dfl 
d/i 



which yields 



1 



5^ — 



~diF 



Sji 



dQg 
d^i 



I 



dpi, 
dQg 
dp 



d^Vtr 



(B23) 



dp^ 



(B24) 



where all quantities are evaluated at the mean field value 
p = p. 

Appendix C: Convergence scheme in the 1/7V 
expansion 

In this appendix we develop a convergence scheme for 
including the effects of 1/A^ corrections to £7. The Gaus- 
sian part of the thermodynamic potential is given by 



^9 = ^E 1^(^11(9)^22(9) - y^lM) 



2/3 



(CI) 



»9i,q 



where 






k 


^iqi - E - E' 


v'^v'^ ' 


iqi+E + E'^ 


Mi2(g) = y \vuv' 

1^ 


( 1 


1 


\iqi+E + E' 


iqi-E- E' 



The expression for Vtg in Eq. (CI) is ill-defined in the 



absence of convergence factors. We remind ourselves here 
that Mil and M22 have different convergence factors. 
We next split ln[MiiM22 - M^j] as 



ln(Mii)e" 



-iqiQ-^ 



ln(M22)e 



-49,0+ 



In 1- 



M?2 



M11M22 

(C2) 
Note, the last term does not need a convergence factor. 
We next note that M22(— (?) = Mii(g) and so upon sum- 
mation over positive and negative values of q, M22 can 
be combined with Mn to give 



rj„ 



^E'" 






/'MiiDetM(g)')e+^^'°" 



(C3) 



We have now endowed U,g with a single convergence fac- 
tor and in order to remove the log divergence at large 
U we now simply add to U,g a term ^- ^ ln([/2)e"'"^'^'° , 
which has no poles or singularities in the left-half plane 
and therefore contributes nothing to the Matsubara sum 
except rendering fig finite in the limit U ^- 00. 

We thus redefine M — > t/M and split the sum in Eq. 



Su ('B22~) (F16 1 into two parts: 



n. 



-y 



In 



Mil 
M25 



„+i9,0+ 



InDetM 



(C4) 



Again, the last term in Eq. (C4) is manifestly convergent 



and hence does not require a convergence factor. How- 
ever, the first term is ultra-violet divergent since in the 
large iqi limit it goes like 



Mil 
M 



1 



2a(q) 



22 



tqi 



(C5) 



where Q;(q) 



c^Ek(' 



2 /2 



V V ). To regulate 



the offending term, we subtract and add a term 



any real number to obtain 



ig+igiO where a is 



-y 

2/3 ^ 
q,«9i 



ln( 



Mil 



a{q) 



1 



1 



M22' '"' \iqi +a iqi - a 

+ ^II"('^) i^T- + -^— 
2(i ^-^ \iqi + a iqi-a 



TC6) 



hi9;0+ 



Now the first term is explicitly convergent and hence the 
convergence factor can be dropped. However, without 
the convergence factor, the first term is an odd function 
of qi and hence the Matsubara sum gives zero. So, we 
are left with only the second term, where the Matsub- 
ara sum can again be converted into an integral along 
the imaginary axis. We can further analytically continue 
on the left half of the complex plane (contribution from 
right half is zero at T = due the Bose distribution func- 
tion) and close the contour counter-clockwise to enclose 
the only singularity at z = —a. The contour integration 
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gives {U/2) X^kl""^"^^) ^'^^ ^hs second term and thus the 
correct form of ilg is given by 



f^„ 



^El-DetM(g) + ^^(^ 



(C7) 



««i,q 



In the next appendix (Appendix pi) we outhne the nu- 
merical steps for evaluating Qg. 



Appendix D: Numerical evaluation of Qg 

The first step in the calculation of the Gaussian part of 
the thermodynamic potential is to solve the gap equation 
for a given chemical potential. Since we do not know the 
analytical form of the number equation once we include 
Gaussian fluctuations, we work in the grand canonical 
ensemble and obtain A(/x) from equations (BIO I. We 



next numeri cally compute rig[/x, A(yu)] using the formula 
in equation (B18). All the 3 momenta sums are over the 
entire Brillouin zone for a 20 x 20 x 20 lattice and have 
an implicit factor of total number of lattice sites in front. 
The Matsubara sum over the imaginary frequencies iqn 
is computed along the imaginary axis for each q mode. 
The integral in equation! B18 1 is split as follows: 



ln(DetM)(iy 



\n(DetM)dy^ 



F{y)dy 

(Dl) 
where the first integral on the left hand side is computed 
numerically and the second integral is evaluated analyt- 
ically using the large y asymptote of the integrand. The 
function F(y) is given by 



^ jfc„,k 



u\'^+v\'^){E + E') 



(D2) 



Here one has to be careful about the integrable log- 
divergence at q = (0, 0,0),y = coming from Gold- 
stone's Theorem. To take this into account we expand 
the integrand for q = (0,0,0) and small y and obtain 



ln(DetM(0,y)) w ln{Ky^), where K 



c? = 2U (l-C^E 



u'^ + v^ 
2E 




u'^ + v^ 



(D3) 



^E 



4 I 4 
-U + V 

4^2 



9 



=-2uy: 



AE^ 



Appendix E: Comparison of our 1/N expansion at 
T = with HQJ (VSR) 

In this appendix we show the equivalence of our 
method for obtaining the 1/N corection to the chemical 
potential at T = with the one obtained by [TUj (VSR). 
We start with equation (B24) and write the numerator 
as 



dQn 



dQo 



dfla dAr 



dfi d^ 9Ao dfj. 



(El) 



We next evaluate the denominator of equation ( B24 ) 

d^n 



d^Qn 



dy? 9/i2 



d^nn 



dAodfiJ dfi 



dAn 



(E2) 



Putting equations (El E2) in equation (B24) we obtain 



J/i 



5/i 



9Af 



dn„\ ( d^no 



aAr 






dAl 



\dAodii 
[dAodfi 



I 
(E3) 



which is the same as equation (3.31) of [TO]. However, it 
should be emphasized that the respective 1/A^ corrections 
to the gap parameter within our theory and within VSR 
are nevertheless different and can be traced back to the 
question of feedback discussed earlier. 



Appendix F: Details of the superfluid density 
calculation for large-TV 

Here we extend our large- A'^ formalism to the calcula- 
tion of the superfluid density P^ in a lattice. We start by 
putting a (time-independent) phase twist Q(x) = Q-x^ = 
di on the order parameter 



A{x) 



'^A{x) 



(Fl) 



b — g and in the expression ( B4 ) for the Nambu-Gorkov propagator. 



which is used to calculate the full action S<s,^a- We shall 
assume that the phase difference is a constant 9 for any 
two neighboring sites along the (arbitrarily chosen) x- 
direction; thus, we have the relation Q = Q^ = 0/a 
where a is the lattice spacing and the superfluid density 
is 



= Ama 






(F2) 



We note that the terms independent of y in the expres- 
sions for a, h and g cancel due to Goldstone's theorem 
and the term linear in y cancel due to symmetry. The 
integrand ln(DetM(0, y)) for q = is then integrated 
between limits and a small value oi y = ys- The rest of 
the integral for q = is evaluated numerically between 
ys and yc, and analytically between limits y^ and cx) using 
the asymptotic form F{y). 



where the mass in the lattice is the combination m — 
l/(2ia2), i.e. the effective mass for fermions at the bot- 
tom of the band. 

We can remove the phase twist from A(a::) by applying 
a unitary transformation to the Grassmann fields of the 
form 



V'ia(T) = XJ^'4>^a{T) 



(F3) 
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with 



_ e-«./2 



(F4) 



so that the inverse Green's function in (B4) now becomes 
Gr^\T, t') -^ U,Gr^\T, t')V]. As can be easily verified, 



this leaves the form of G^ {t,t) unchanged from (B4) 
except for the hopping term, which gains the phase dit- 
ference t ^ t exp{zLi{0j — di)/2) in the first (second) 
diagonal element. Transforming to the reciprocal lattice, 
we obtain an inverse Green's function of the form 



Go^ = 



(ik„ 



Ck+Q/2)e" 
Ao 



-ifejO^ 



Ao 
(ifc„ + Ck-Q/2)e 



+ikiO+ I 

(F5) 

Taking the limit of small 9, we see that this corresponds 
to shifting the Matsubara frequencies and the energy dis- 
persion respectively as 



l/iXi'ry 1/ fXi'y-i 



a = 6 



{de^/dk^) 



2a 
8a2 dkl 



(F6) 



The effective action at a fixed 9, from which the saddle 
point condition is derived, satisfies 



5o(Ao;^,g) ^ Ag 1 
U p 



NP 



--^TrlnG„i(fc) 

ikji ,k 



(F7) 



iV/3 



2 2^^ Ejdkl 



where we have expanded the energy dispersions to 
quadratic order in 9 and used the fact that the shift in 
the Matsubara frequencies is not important once they 
are summed over as long as 9 is small. The saddle point 
condition 5Sq/51S.q = yields the small 9 expansion 
Ao(m, d) = Ao(m, 0) + a{^l)9'^, where 



a(/i) 



3 T.^{^^^,./^kl)i^^|El 



(F8) 



8Aoa2 Ek(l/i5^) 

and all quantities are evaluated at = 0. 

1. Mean Field Superfluid Density 

The mean field thermodynamic potential per fiavor at 
T = for a system with a phase twist is given by 

1 



n^{^J^,9) = —So{/^o{^i,9)■^i,9) 



(F9) 



Using that the Aq dependence on 9 is obtained from the 
saddle point condition (55o/(5Ao = we can obtain the 
superfiuid density as 

1 fd'^S^" 



NP \ 89^ 



/J, A 



Eh 



Ev 



cos{kxa) 



(FIO) 
(Fll) 



2. Calculation of n^ including Gaussian fluctuations 

In order to include Gaussian fluctuations, we need to 
calculate the Gaussian part of the action Sg{Ao;9,^) in 
the presence of the phase twist 9. The inclusion of the 
effects of Gaussian fluctuations in the calculation of the 
superfiuid density follows the same methodology used in 
section [VC| The thermodynamic potential per flavor to 
first order in 1/A^ is of the form 



n(M,0) 



^oil^,9) 



1 

'np 



ng{f,,9). 



(F12) 



where again Qg{fi, 9) = Sg{/^o{n, 9);^, 9)/p. 

Before we give an explicit expansion of the superfiuid 
density in orders of 1/A^ on a lattice, let us go back to 
the continuum limit and outline the calculation of the 
1/N corrections to rig in the continuum. This would 
hopefully elucidate some of the technical points differ- 
entiating a lattice calculation from the continuum. To 
this effect we prove that for a translationally invariant 
system, the relation Ug — n, is respected even at the 1/N 
level. For an energy dispersion Ck = k2/2?7i, the shift in 
single particle energies due to the introduced twist (see 
second line in Eq. F7 ) is only a constant and can be in- 



corporated into the chemical potential. Since, the phase 
twist is uniformly distributed across the system the or- 
der parameter transforms as: A(x) — ?> A(2:)e*'^'' and the 
Green's function transforms as 



G,y (ifc„, k, /i) -^ Gij(zfc„, k, n-Q 



im) 



(F13) 



«'7/,k 



where the ifc„ = ifc„ — k.Q/2TO are the Doppler shifted 
Matsubara frequencies. Then, 'M.ii{iqi,q;Q) ~ 1/U + 

.G,,{^k„ - ^,k;M - E)Gii(^fcn - ^"^^^ + 

q;M-£) = Mn(*g,-'i^,q;M-£). The 
effects of the phase twist at T = is therefore to shift 
the contour of integration for the Matsubara sum by an 
amount proportional to Q along the real axis and to shift 
the chemical potential /z by a constant amount Ci^/Sm. 
In the limit Q — >■ 0, the shift of the contour of integra- 
tion to the right keeps the Matsubara sum invariant and 
hence the phase twist enters the thermodynamic poten- 
tial only through a shift of the chemical potential. This 
means that the saddle point condition in presence of the 
phase twist remains unchanged from the one in absence 
of the same. Further, fl{Q) only contains terms in powers 
of Q^ and therefore we obtain. 



Am 



&H}g 

5Q2 



(Klg 



(F14) 

By the same logic, Am{dfto/dQ'^)jz^Q^o = -(9r2o/5/i)jr. 
Since, the number equation is given by d(i2o + 
Q,g/N)/dii — —n we obtain 



rig = Am 



dno 



1 

N 



dn. 



dfi 



n 



(F15) 
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We next consider the lattice and write the 1/N correc- tions to the thermodynamic potential as follows 

+ ^ng{p,^^{^i,e),e) (fi6) 

Further, fi ^ JI + 6fi/N, which then combined with Eq. 
( F2 ) yields the following expansion for the superfluid den- 



sity 





-[d0^ ) 


1 d fd^no\ 1 fd^ng\ 


>^. 


' d fd^flo^ 
dfi \ de^ , 





(F17) 



where d/dfi — d/dfi + (9/x/(9Ao)9/5Ao. Note the pres- 
ence of an explicit d/dAo derivative which was absent in 
the expression for n^ because of the saddle point condi- 
tion 9r2o/9Ao = 0. Eq. |F17 gives the expression for n^ 
that was used in section fVIl 



Appendix G: Critical temperature using Hartree 
shifted NSR 

In this appendix we shall use the Nozicrcs and Schmitt- 
Rink approach [3] for calculating the critical temper- 
ature, with the modification that the single particle 
Green's function Gq now includes a Hartree shift which 
we call S. We shall work in the grand canonical ensem- 
ble at a fixed fj, and hence E is a function of n and T. 
We approach the transition from above Tc and look for 
the divergence of the t-matrix. This gives us a relation 
between Tc and /i 



i=5:^tanh( 



/^^ 



U 



) 



(Gl) 



Here the Hartree shift is contained in fk = ek — M + ^■ 
Since we are working in the grand canonical ensemble the 
filling fraction would depend on the value of n we choose. 
At a fixed temperature this dependence is given through 
the number equation 



^niT,^l) 



d^i 



(G2) 



The Hartree shift, which depends on the filling fraction, 
is then given by 



J:{^^,T) = -n{^l,T)U/2 



(G3) 



In order to implement the number equation ( G2 ) wc pro- 



ceed as follows: For a given U and /i, we calculate Tc and 
T,{fi,Tc) by simultaneously solving Eqs. (Gil and (G3|. 



I 

With these values of Tc and E, we evaluate il(/x, T = Tc). 
Next, keeping the temperature fixed we change fi to ii+Sfi 
and evaluate J:{fi + S^i,T = Tc) from Eq. ([^. This lets 
us evaluate fl{fj, + 6fi, T = Tc). The number equation can 
then be written in the form 



n{^i + Sfi,T = Tc) - n{n, T = Tc) 
6^1 



(G4) 



We next give an explicit formula for the thermody- 
namic potential 17 = 51o + ^g- In presence of the Hartree 
shift, the mean field thermodynamic potential fig is given 

by 



no 



EMI 



-/3?k 



E2 
U 



(G5) 



As a check, note that setting O = Oq in Eq. (G2) the 



above form for Oq gives us the familiar mean field number 
equation 



E 



exp(/?^k) + 1 



(G6) 



To obtain the thermodynamic potential upto Gaussian 
order for T > Tc, we note that Uk = 1 and Vk = and 
hence M12 = M21 = 0. Therefore 



"ff = ^EM(Mii(<z))-(Mn(<z) 



/3 



1)] 



(G7) 



q,*9i 



where, we are justified to drop the convergence factor 
in the last line of Eq. (G7| since (Mn — 1) takes out 
the leading l/iqi piece from In(Mii) (see Appendix [A| 
for more details). For T > Tc, Mn(g) = 1 + C/^k(l - 
/ ~ f')/i^Qi ~ "? ~ 0, where iqi — i2-Kl/ ji are the Bose 
Matsubara frequencies and f{£) — l/[exp(/3^) -I- 1] is the 
Fermi distribution function. We numerically calculate ilg 
using Eq. ( G7 ) and set up the number equation using 
the procedure outlined above. In order to obtain the n 
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the procedure is repeated for various 

We 



dependence of T^ 

values of /i. The results are plotted in Figure ( 12 ) 



notice that there are quantitative differences between the 
HNSR and large--/V results. At this stage, we do not 
understand why the Tc from HNSR is lower than the Tg 
from large- A^ theory. 



FIG. 12. Figure shows a comparison of the critical temper- 
atures obtained within Hartree shifted NSR and large A'' re- 
spectively. The filling fraction n — 0.5. The maximum Tc 
for the Hartree shifted NSR and the large A'^ theory can be 
compared to the results of Tamaki et al. (T™^'' ~ 0.66i for 

n = o.5) Eg 
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